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ABSTRACT 

In this paper we present the results of time-dependent simulations of dipolar ax- 
isymmetric magnetospheres of neutron stars carried out both within the framework 
of relativistic magnetohydrodynamics and within the framework of resistive force-free 
electrodynamics. The results of force- free simulations reveal the inability of our numer- 
ical method to accommodate the equatorial current sheets of pulsar magnetospheres 
and raise a question mark about the robustness of this approach. On the other hand, 
the MHD approach allows to make a significant progress. We start with a nonrotat- 
ing magnetically dominated dipolar magnetospheres and follow its evolution as the 
stellar rotation is switched on. Wc find that the time-dependent solution gradually 
approaches the steady state that is very close to the stationary solution of the Pulsar 
Equation found by Contopoulos et al.(1999). This result suggests that other station- 
ary solutions that have the y-point located well inside the light cylinder are unstable. 
The role of the particle inertia and pressure on the structure and dynamics of MHD 
magnetospheres is studied in details, as well as the potential implications of the dis- 
sipative processes in the equatorial current sheet. We argue that pulsars may have 
differentially rotating magnetospheres which develop noticeable structural oscillations 
and that this may help to explain the nature of the sub-pulse phenomena. 

Key words: stars:pulsars:general - stars:winds and outflows - relativity - MHD 



1 INTRODUCTION 

Pulsar magnetospheres are complicated electrodynamic sys- 
tems that involve the acceleration of charged particles by the 
rotationally induced electric field, the electron-positron pair 
production, and the bulk outflow of magnetospheric plasma 
through the light cylinder. The pair production presumably 
occurs on the scale which is significantly smaller compared 
to the magnetospheric scale (the last is determined by the 
light cylinder radius vji c = c/fi). This suggests that that 
the global structure of pulsar magnetospheres may not de- 
pend on the details of the pair production mechanism which 
may enter the problem only in the form of boundary condi- 
tions. The simplest assumption is that the plasma supply is 
sufficiently plentiful to allow the MHD description through- 
out the whole magnetosphere. Although this may be not so 
in the case of slowly rotating or weakly magnetized pulsars 
(Hibschman & Arons 20011, it still makes sense to regard 
the MHD solution as a zero approximation model and a 
good reference point. 

Although pulsar magnetospheres are generically three 
dimensional it is still reasonable to start with a much simpler 



axisymmetric case as it involves most of the basic physical 
elements of the general case. Michelfl9731 and Scharlemann 
& Wagoner 1119731 first analysed the structure of the sta- 
tionary axisymmetric pulsar magnetosphere in the limit of 
force-free electrodynamics where both the pressure and iner- 
tia of magnetospheric plasma are assumed to be vanishingly 
small. They found that in this case the problem is reduced 
to the following second order PDE that is now called the 
Pulsar Equation, 



(1 



dH_ 1 <9* d 2 t> 

dx 2 x dx dy 2 
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where x — uj/vji c and y — z/vji c are the cylindrical coordi- 
nates normalised to the radius of the light cylinder and $ 
is the magnetic flux function. This equation also involves 
the poloidal current function, A(9), that is unknown. If 
this function is somehow specified the pulsar equation be- 
come a linear elliptic equation. However, the question of 
what actually determines the electric current function is not 
that obvious and is still debated ( Mcstcl 1999 Bcskin 2005). 
The physical meaning of these functions is quite simple: 
* = $/27r, where &(x,y) is the magnetic flux through the 
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axisymmetric circular loop given by x =const, y =const; 
A — 21 /c where I(x,y) is the total electric current flowing 
through the same loop. It is also worth mentioning that if we 
set A<f,(x — 0) = then $ = A^, and that in a steady-state 
A = H^,, where A^ and H$ are the covariant azimuthal com- 
ponents of the vector-potential A and the magnetic field H 
in the non- normalised basis of cylindrical or spherical coor- 
dinates. 

Another important feature of this equation is the math- 
ematical singularity at the light cylinder. It was argued 
by Ingraham ll973t that the condition of smooth passage 
through this surface together with the appropriate bound- 
ary conditions determine the unique electric current function 
of the pulsar equation (This property makes the problem 
somewhat similar to the classical eigenvalue problem in the 
theory of differential equations.) Ingraham|1973|l even pro- 
posed an iterative algorithm for finding this function. In the 
same year Michel(iQ73) did actually find an exact solution 
to the pulsar equation in the case of split-monopole mag- 
netic field that passed through the singular surface both 
continuously and smoothly. Although Michel did not use 
Ingraham's approach his result supported Ingraham's idea 
and raised expectations for finding smooth continuous so- 
lutions in more complicated realistic cases, like for example 
the magnetosphere of a rotating dipole, as well. (The recent 
observations of magnetospheric eclipses in the binary pul- 
sar J0737-3039 have confirmed that the dipolar model may 
indeed be quite accurate jLyutikov fc Thom pson 2005J.) 

The electric current of Michel's solution has the same 
sign within the magnetosphere with the exception of the 
equatorial current sheet that provides current closure. Fur- 
ther analysis by Ingraham <ll973H and Michel <ll974H sug- 
gested that if the force-free solution for an axisymmetric 
rotating dipole existed then far beyond the light cylinder it 
should still closely resemble the split-monopole one. How- 
ever, in the near zone the structure of the dipole magneto- 
sphere was expected to be qualitatively different. Indeed, 
close to the star the effects of its rotation are relatively 
small and so must be the difference between the solutions 
for the rotating and the non-rotating magnetospheres. As 
the result, the so-called "dead zone" should exist where the 
magnetic field lines remain closed and the magnetospheric 
plasma co-rotates with the star. 

However, the problem of finding the self-consistent 
global solution for dipolar magnetospheres turned out to be 
rather involved, even for the axisymmetric case of aligned ro- 
tator. For a very long time, the only attempts to construct 
such solutions involved utilisation of prescribed electric cur- 
rent functions. Such solutions would either exhibit kinks at 
the light cylinder or violate the conditions of force-free ap- 
proximation at some relatively short distance from it, e.g. 
JMichel 19821 IBeskin et al. 19931 . The break through came 
only recently then Contopoulos et al. (1999, hereafter CKF) 
have finally managed to solve the problem numerically. In 
this solution the dead zone continues all the way to the light 
cylinder where the so-called "Y-point" appears in the mag- 
netic field structure whereas beyond the light cylinder it has 
the same topology as the split-monopole solution including 
an infinitely thin equatorial current sheet. The return cur- 
rent of the equatorial current sheet splits at the Y-point 
into two currents sheets that flow along the surface separat- 
ing the dead zone from the open field zone. An additional 



finite width layer of return current exists just outside of the 
dead zone and around the equatorial plane. 

Uzdensky { 2003 ) carried out asymptotic analysis of the 
force- free solutions in vicinity of the Y-point. It has turned 
out that in the presents of separating current sheets, sim- 
ilar to those found in ( |Uontopoulos et al. 1999| , the elec- 
tromagnetic field becomes infinitely strong at the Y-point. 
He argued that this is unphysical and questioned correct- 
ness of the CKF-solution. As an alternative Uzdensky ( 2003 1 
considered the case with no current sheets, that is the case 
where the electric current returns within the equatorial layer 
of the open field zone only. In this case the electromagnetic 
field does not diverge at the Y-point but just outside of this 
point the electric field becomes stronger than the magnetic 
field signaling a breakdown of the force-free approximation. 
No such complications were found for the Y-point located 
inside the light cylinder. In fact, in this case the separating 
surface crosses the equator at the right angle so that the 
Y-point turns into a "T-point". 

Gruzinov( 2005 1 also analysed the force-free solution in 
the vicinity of the Y-point located at the light cylinder in 
the presence of the current sheets and confirmed that in this 
case the electromagnetic field diverges at the Y-point. His 
analytical solution exhibits the angle of 77°. 3 between the 
equatorial plane and the separating surface at the Y-point. 
In order to verify this he repeated the CKF-calculations with 
much higher spatial resolution and the results seemed to 
confirm the development of a singularity with the expected 
inclination angle of separatrices. 

The next potentially very important step was made by 
Goodwin et al, 12UU4(l who realised that the dead zone does 
not have to extend all the way to the light cylinder but can 
be much smaller. Such solutions select their own electric cur- 
rent functions, smaller dead zones corresponding to weaker 
currents. In fact, the first solutions of the pulsar equation 
which described both the equatorial current sheet and the 
dead zones that were located deep inside the light cylinder 
were found by Lyubarskii( 1990). However, these solutions 
utilized a prescribed electric current function and, as the 
result, exhibited the breakdown of the force-free approxi- 
mation - at a certain distance from the light cylinder the 
electric field became stronger than the magnetic field. Good- 
win et al.(2004| also included finite gas pressure inside the 
dead zone and showed that this allows solutions that remain 
nonsingular at the Y-point even when this point is located 
on the light cylinder. In this case their solution gives the 
inclination angle of the Y-point of 56°. 5. However, the in- 
ertia associated with the gas pressure has not been taken 
into account and thus the self-consistency of such approach 
is questionable. 

This idea of Goodwin et al.(2004 ) was further explored 
by Timokhin ll2U05t who has constructed numerical solu- 
tions for force-free magnetospheres with the Y-point located 
within the light cylinder (with no finite gas pressure in the 
dead zone). He also pointed out that the pulsar spin-down 
rate depended not only on the size the light cylinder but 
also on the size of the dead zone and if the ratio of these 
scales were to evolve with the pulsar age then this could ex- 
plain why the observed braking index of pulsars is smaller 
than the one derived from the models with the dead zone 
extending all the way up to the light cylinder. 

Another interesting twist has been added by Contopou- 
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los(2005Jl, who argued that during the pulsar evolution the 
potential gap separating the star surface from the polar cap 
magnetosphere grows in magnitude leading to a slower ro- 
tation of the open field lines compared to the star and the 
dead zone. In such a case, there is no a single light cylinder 
for the whole magnetosphere. Instead, the dead zone has its 
own light cylinder that is located inside the light cylinder 
of the polar cap. To illustrate this point Contopoulos(2005) 
constructed a set of numerical models of pulsar magneto- 
spheres of this kind. The evolution of the angular velocity 
of the polar cap magnetosphere with pulsar age may also 
result in the lower values of the pulsar's braking index. 

In spite of this remarkable progress the method of pul- 
sar equation has its obvious limitations. It does not ad- 
dress the question of stability of the steady-state solu- 
tions and it does not allow to study possible non-stationary 
phenomena in pulsar magnetospheres. Moreover, the nu- 
merical techniques that have been used to solve the pul- 
sar equation may turn up to be rather difficult to adapt 
to the fully 3-dimensional problem of the oblique rotator. 
The obvious way of approaching these problem is to relax 
the stationarity condition and to solve the original system 
of time-dependent equations. Although the force-free elec- 
trodynamics was used to model various astrophysical sys- 
tems since 1970s the focus was entirely on the steady-state 
equations. Only recently, the time-dependent equations 
have been subjected to a systematic study (lUchida 19971 
IGruzinov 19991 IKomissarov 2002al |Punsly 2001) . As the re- 
sult it has been found that they form a simple hyper- 
bolic system of conservations laws in many respects sim- 
ilar to relativistic MHD but only with the fast and the 
Alfven hyperbolic waves (IKomissarov 2002at . Thus, a va- 
riety of standard numerical methods can be used to deal 
with these equations. The very first numerical simulations 
of this kind, of the monopole magnetospheres of black 
holes, seemed to confirm the suitability of this approach 
(Komissarov 20011. However, further applications to some- 
what more complex magnetic configurations revealed its 
limitations as the force-free approximation would often 
break down following development of strong current sheets 
omissarov 2002a] |K omissarov 2002b] ISpitkovsky 2004| 
lAsano et al. 1999ll . These results suggested to look for a 
more general mathematical framework that would allow to 
handle current sheets via introducing new channels of back 
reaction of plasma on the electromagnetic field. One of the 
most suitable options is the framework of resistive relativis- 
tic MHD. However, this approach remains rather poorly 
studied even now. Ideal relativistic MHD is a somewhat eas- 
ier option and thanks to the recent effort by several groups, a 
significant expertise has been obtained in developing numer- 
ical methods for this system. This approach allows to take 
into account both the thermodynamic pressure of plasma 
heated in the current sheets but the dissipation of electro- 
magnetic energy is entirely due to numerical resistivity that 
is not fully satisfactory. 

Another alternative is resistive force-free electrodynam- 
ics with physical or artificial resistivity. If the current sheets 
of pulsar magnetospheres are indeed dissipationless then this 
approach might work provided that the utilized Ohm law al- 
lows evolution towards a dissipationless force-free state. The 
recent time-dependent simulations by A.Spitkovsky show 
that this approach may indeed be quite productive (the 



results were presented last summer at the conference on 
"Physics of Astrophysical Outflows and Accretion Disks", 
KITP, Santa-Barbara, 2005.) 

In this paper we report the results of new time- 
dependent axisymmetric simulations of rotating dipole mag- 
netospheres of neutron stars within the frameworks of resis- 
tive force-free electrodynamics and ideal relativistic MHD. 
We only consider the case where the whole magnetosphere 
rotates with the same angular velocity and thus our re- 
sults are not relevant to the model of Contopoulos(2005 ). 
This model will be a subject of a separate study. Note, that 
throughout the paper the graphic data are shown in the di- 
mensionless units where the magnetic dipole moment fj, = 1, 
the angular velocity of the star SI, and the speed of light 
c = 1. Hence the cylindrical radius of the light cylinder 

UJl c — 1. 



2 COMMON FEATURES OF THE 
SIMULATIONS 

Using the standard notation of the 3+1 approach the metric 
form of a general space-time can be written as 

ds 2 = (P 2 ~ ct 2 )dt 2 + 2Pidx i dt + " 1i] dx i dx i , (2) 

where 7^ is the metric tensor of "the absolute space", a 
is the lapse function, and (3 is the shift vector. For most 
purposes in physics of pulsar magnetospheres the flat space 
approximation suffices and one can enjoy the benefits of a 
global inertial frame where a = c and f3i = (see however 
Beskin 1990 and Muslimov & Tsygan 1990). This is exactly 
what was adopted in the numerical models of pulsar mag- 
netospheres described in the Introduction and in this study 
as well. However, the computer codes used in the simula- 
tions are designed to work with a rather general axisym- 
metric and stationary space-times. Moreover, in the case of 
oblique rotator the frame rotating with the star seem to be 
more suitable as the solution may become stationary in this 
frame. In such a case (3 is no longer vanishing and in the 
basis of spherical spacial coordinates 

Pt =08=0, fa = cfisin 2 0r 2 , 

and hence 

2 2 i ry2, 2 r~\ 2 ■ 2 r\ 2 

a = c , and p = c il sin or , 

where SI is the angular velocity of the frame. 

In our simulations, the computational grid covers the 
axisymmetric domain (r,6) = [ri n ,r ou t] x [0,7r] so no sym- 
metry condition is enforced at the equatorial plane. The cell 
size in the r-direction is such that the corresponding phys- 
ical lengths in both directions are equal, Ar = rAO. The 
"radiative" outer boundary, r = r out , is always located so 
far away from the star that the light signal does not cross 
the computational domain by the end of the simulations - 
this ensures that waves produced near the star do not get 
reflected of the outer boundary at any rate. 

As the Courant stability condition requires At < Ar/c, 
the suitable time step for the outer part of the computa- 
tional domain is much larger than the one for its inner part. 
This allows us to reduce the computational cost of simula- 
tions via splitting the computational domain into a set of 
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rings such that the outer radius of each ring is twice its in- 
ner radius and advance the solution for each fcth ring with 
it its own time step, Atk, such that Atk+i = 2 Atk- Thus, 
one integration step of ring k corresponds to two integra- 
tion steps of ring k — 1, four integration steps of ring k — 2 
and so on. As the result, the outer regions of the compu- 
tational domain are progressively less expensive in terms of 
CPU time. This approach has already been successfully ap- 
plied in recent MHD simulations of Pulsar Wind Nebulae 
l |Komissarov fc Lyubarsky 2004b] > . 

The initial solution describes a nonrotating magneto- 
sphere with exactly dipolar magnetic field, 

= 0; 

B f = 2p cos 9/r 3 ; (3) 
B e = fi sin 9/r 3 , 

where fj, is the magnetic dipole moment. At time t — 
the stellar rotation is suddenly switched on via an appro- 
priate inner boundary condition. The simulations are then 
proceeded till it becomes clear whether the time-dependent 
solution relaxes to a steady state on scales comparable with 
the light cylinder radius or not. Steady-state solutions that 
are found in such a way are automatically stable to axisym- 
metric perturbations with wave-lengths exceeding the cell 
size. 
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Figure 1. Representative force-free solution. The contours show 
the magnetic flux surfaces and the colour image shows H j, . Notice 
that the magnetic field lines are closed even beyond the light 
cylinder, ro; s = 1. 



3 ELECTRODYNAMIC MODEL 
3.1 Equations 

Following Komissarov( 2004a ) we write vacuum Maxwell's 
equations as 

d i (^B i ) = 0, (4) 

(l/c)d t B i + e» k d j E k =0, (5) 

driVlD 1 ) = A*k^k, (6) 

- (l/c)d t D* + , ■'.t,!!, = (4tt/c)J\ (7) 

Here 7 = det(7ij) is the determinant of the metric tensor of 
the absolute space, e-ijh = \/7 e »J* * s ^ ne Levi-Civita tensor 
of the absolute space (ei23 = 1 for right-handed systems and 
£i23 = — 1 for left-handed ones.) The electric field, E, and 
the magnetic field, B, are defined via 



01 

Ei — ~€-i-jk E 

2 

and 



■jk 



(8) 



(9) 



where *E y,v is the Faraday tensor of the electromagnetic 
field, which is simply dual to the Maxwell tensor _F M ": 



2 



where 



(10) 



(ii) 



is the Levi-Civita alternating tensor of spacetime. 

The differential equations I4I7H are supplemented with 
the following constitutive equations (Komissarov 2004at : 



cE u = aDk + ekijpBi, 
cH k = aBk - ekij^D 1 . 
J k = (a/c)j k - K,p k . 



(12) 
(13) 
(14) 



Note, that 1) the components of all vectors and tensors 
appearing in these equations are measured in the non- 
normalised coordinate basis, {di}, of spacial coordinates; 2) 
D, B, j, and re are the electric field, the magnetic field, the 
electric current density, and the electric charge density as 
measured by a local observer at rest in the absolute space. 
The 4- velocity of this observer is 

n = ~(c,P'); 
a 

In the inertial frames of flat spacetime E = D and B — H. 

The final equation that is needed to close the system is 
the Ohm law. In strongly magnetized plasma the conductiv- 
ity is no longer isotropic and under rather general conditions 



3 = a \\ D \\ + cr ±- D ± +3d, 



where 



3d = KC " 



DxB 



B 2 



(15) 



(16) 



is the drift current. Unless the magnetosphere has scarce 
supply of electrically charged particles /'charge starved", 
the parallel conductivity is very larger that leads to small 
residual parallel component of electric field, 



0. 



(17) 



However, the strong magnetic field of pulsar magnetospheres 
effectively suppresses conductivity across the magnetic field 
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lines (unless the electric field is even stronger than the mag- 
netic one; this may be the case inside some current sheets.) 
So we can put 



a ± =0 if B 2 >D 2 . 



(18) 

As shown in (JKomissarov 2004a), in this limit we approach 
the approximation of force-free electrodynamics. 

Within current sheets this simple prescription is un- 
likely to hold. On one hand, in the area of high current 
density one may expect strong anomalous resistivity leading 
to significantly reduced cry. As the result, the current sheet 
may become unstable to the tearing mode of reconnection 
<Lyutikov 2003) . Since our intention at this point is merely 
to see if we can construct idealized numerical solutions that 
are more-or-less force-free and thus can be compared with 
the steady state solutions of the Pulsar Equation we will 
ignore this physical effect for the time being. On the other 
hand, the mere symmetry of the aligned rotator magneto- 
sphere suggests that sooner or later the magnetic field of 
force-free solution will become very small in the equatorial 
plane beyond the light cylinder, thus leading to unavoidable 
breakdown of the D 2 < B 2 condition of the force-free ap- 
proximation. In such a case of relatively weak magnetic field 
the conductivity is expected to become much less anisotropic 
and we will assume that 



(7_L = cry > 1 if B 2 < D 2 . 



(19) 

Notice, that the expression Ijlfip for the drift current should 
also be modified in this case as it implies that the drift speed 

DxB 



Vd = c- 



B 1 



(20) 



becomes higher then the speed of light. We have tried several 
modifications for the Vd including 

DxB 

Vd = C m»x{B*,D*) (21) 
which is likely to underestimate va in the current sheet and 

Vd = Cjr^. — if B < D 



\DxB\ 



(22) 



which certainly overestimates it. The actual value of conduc- 
tivity in the current current sheet remain a free parameter. 

To find numerical solutions of Maxwell's equation, in 
particular the solutions that are presented in the next sec- 
tion, we used the Godunov-type numerical scheme described 
in Komissarov( 2004a ). 

3.2 Numerical simulations 

The most striking and at first somewhat perplexing result of 
our electrodynamic simulations of pulsar magnetospheres is 
demonstrated in figure 0- contrary to what is found in the 
steady-state solutions of the pulsar equation the magnetic 
field lines remain close even outside of the light cylinder. 
Moreover, the solution seems to be rather insensitive to the 
details of the model for a± and remains qualitatively the 
same even in the case <jj_ = (In fact, the solutions exhibited 
convergence in the limit a± — > 0.) 

The question of whether the field lines extending be- 
yond the light cylinder should open up or not is in fact rather 
involved. One of the arguments often put forward in its dis- 
cussion concerns the requirement for the speed of charged 



particles that fill the magnetosphere to remain smaller than 
the speed of light. In the drift approximation the charged 
particles move along the rotating magnetic field lines like 
beads on wire - their motion is a composition of the ro- 
tational motion of the field line (the "wire") and the slid- 
ing motion of the particle (the "bead") along the field line. 
Within the light cylinder the speed of the rotational mo- 
tion is less than the speed of light and the bead does not 
have to slide along the wire, but this is no longer the case 
beyond the light cylinder. Here the total speed of the bead 
can remain less than c only if it also slides along the wire 
and only if the wire is twisted in the azimuthal direction (In 
such a case the sliding motion may help to reduce the az- 
imuthal component of the bead velocity.) These conditions 
may well be satisfied everywhere along the open field lines 
but not along the closed ones. Indeed, such field lines can- 
not have an azimuthal component in the equatorial plane 
because of the symmetries of the problem. Thus, one may 
expect spinning up of the charged particles till their inertia 
becomes important and the centrifugal force opens up the 
closed field lines. 

On the other hand, the drift approximation itself may 
breakdown and the charges may start moving across the 
magnetic field well before their inertia becomes dynam- 
ical significant. One can look at this breakdown of the 
drift approximation from another prospective. Beyond the 
light cylinder the electric field of a force-free solution be- 
comes stronger than the poloidal component of the mag- 
netic field and thus stronger than the total magnetic field 
in regions where the azimuthal component of magnetic field 
vanishes. Such a strong electric field is capable of "tear- 
ing" charged particle off the magnetic field lines and driving 
strong electric current across the magnetic field. This is a 
micro-physical argument but there is an equally compelling 
macroscopic one. Full opening of magnetic field lines implies 
an infinitely thin equatorial current sheet with an infinitely 
high electric current density. Even very small but finite re- 
sistivity will destroy this ideal configuration and result in a 
steady-state current layer of finite thickness and closure of 
some of the magnetic field lines in this layer. Thus, the only 
meaningful question is how many field lines will be closed in 
this layer. 

In order to explain the remarkable indifference of our 
solutions to the the value of a± it is instructive to consider 
a much simpler problem of a one dimensional current sheet. 
Here we adopt an inertial frame with Cartesian coordinates 
so that (D = E and H = B). Assuming that B = (0, B v , 0) 
and E — (0,0, E z ) and using the following model for the 
cross-field conductivity 



(To 





if 
if 



E 2 > B 
E 2 < B 2 



(23) 



we find the following solution 




E" 



X < —l/<7o 
— 1/(70 < X < 1/(70 
X > 1/(70 

B. 



This solution exhibits a number of interesting features. First 
of all it is stationary. Secondly, the magnitude of oo effects 
only the width of the current sheet. This may explain the 
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observed convergence of the magnetospheric solutions in the 
limit <tj_ — * 0. Finally, the electromagnetic energy flows into 
the current sheet with the speed of light and disappears in- 
side of it. This property clearly exposes the nature of Ohm's 
law 1151 with prescription 123H or similar and outlines the 
limitations of its applicability. Indeed, in real plasma the 
Ohmic dissipation would cause strong heating and increase 
of the gas pressure in the current sheet. This pressure would 
slow down the inflow of plasma into the current sheet and 
significantly modify its structure. However, this factor, as 
well as the plasma inertia, are completely ignored in our 
version of Ohm's law. 

The fact that our numerical magnetospheric solution 
remains qualitatively the same even in the limit <t_l = 
suggests that the numerical resistivity of our scheme takes 
over in this limit. Different numerical schemes have different 
dissipative properties and this may explain why the recent 
electrodynamic simulations of Spitkovsky (which have been 
presented at various astrophysical meetings) exhibit opening 
up of the field lines beyond the light cylinder and develop- 
ment of the dissipationless equatorial current sheet. 

The properties of the current sheets of real pulsars 
are determined by a number of competing micro- physical 
processes that might be rather difficult to account for in 
a macroscopic model. Our experiments with the simplistic 
generalised Ohm law H15H22H show that it facilitates devel- 
opment of strongly dissipative current sheets. Since the dissi- 
pated energy of the electromagnetic field is not stored in any 
dynamical component and simply vanishes from the system 
this model could only be relevant for magnetospheres with 
effective radiative cooling. The strong gamma-ray emission 
from some pulsars may be interpreted as an indication of 
strong radiative current sheets (see the Discussion). Unfor- 
tunately, such emission is not a common feature of pulsars. 



4 MHD MODEL 

The difficulties that we have encountered in dealing with 
current sheets in force-free electrodynamics have forced us 
to return to the framework of full relativistic MHD even if 
the pulsar magnetospheres are expected to be magnetically 
dominated everywhere else and in spite of the fact that the 
system of MHD equations becomes stiff in this regime. The 
MHD approximation takes into account not only the gas 
pressure but also its inertia which may become important 
both far away from the star, in the wind zone, and near the 
light cylinder. Ideally we should also incorporate a physical 
model for electric resistivity which would allow us to con- 
duct a proper study of the dissipation within the equatorial 
current sheet but it makes perfect sense to start with the 
simpler framework of ideal MHD. 

4.1 Equations 

The system of ideal relativistic MHD includes the continuity 
equation 

dtiot^jpu 1 ) + d^a^pu 1 ) = 0, (24) 

where p is the rest mass density of matter and u u is its 
4- velocity; the energy-momentum equations 



d t {a^T\) + diia^Ti) = ~d„(g a/3 )T ae a^j, (25) 

where is the total stress-energy-momentum tensor; the 
induction equation, 

{l/c)d t (B l ) + e iik d j (E k ) = 0, (26) 
and the divergence free condition 

8 i (^B i )=0. (27) 

The total stress-energy-momentum tensor, T^" ', is a sum of 
the stress-energy momentum tensor of matter 

T^ ) =wu»u»-pg>"', (28) 

where p is the thermodynamic pressure and w is the enthalpy 
per unit volume, and the stress-energy momentum tensor of 
the electromagnetic field 

T% = ^ {P^F% - \{F^F aP )g^) , (29) 

where F vtl is the Maxwell tensor of the electromagnetic field. 
In the limit of ideal MHD 

Ei = e ijk v i B k /c, (30) 
where v* = u l / u * i s the usual 3- velocity of plasma. 

4.2 Numerical method 

The MHD simulations were carried out using a Godunov- 
type scheme that is described in Komissarov( 1999 2004b). 
However, we had to introduce a number of additional fea- 
tures in order to overcome a number of challenging problems 
specific to the case of highly-magnetized plasma. 

The MHD equations become stiff in magnetically dom- 
inated domains and this is exactly the case for the main 
volume of pulsar magnetospheres where the energy density 
of matter is many orders of magnitude less than the one of 
the electromagnetic field. There is no hope of reaching such 
conditions with our numerical method. However, the elec- 
tromagnetic part of the MHD solution should not be very 
different even when the energy ratio is artificially increased 
up to 0.1-0.01 - errors of order of few percent are quite ac- 
ceptable at this stage of investigation. In fact, the results 
of force-free and MHD simulations of the monopole mag- 
netospheres of black holes strongly support this conclusion 
( Komissarov 2001 Komissarov 2004b I . In those MHD simu- 
lations, additional plasma was pumped in the regions where 
its energy density had reached a certain lower limit. Here, 
we apply a similar trick. The actual condition is 

wW 2 > a {1) B 2 , (31) 

where W is the Lorentz factor of plasma and am is a small 
constant (we used am = 0.01.) When this condition is bro- 
ken we increase both p and p in the same proportion so that 
wW 2 = a (1) B 2 . 

However, the dipolar magnetospheres are more chal- 
lenging compared to the monopole ones because of the faster 
decline of magnetic field strength with the distance from 
the star and the existence of dead zones. Within the dead 
zones the magnetospheric plasma is supposed to be in static 
equilibrium. Thus, the component of the centrifugal force 
acting along the magnetic field lines has to be balanced by 
some other force. It has been argued that the dead zone 
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Figure 2. The evolution of spinned up dipole. This figure shows the solution at t = 5 (left panel) and at t = 20 (right panel). The 
contours show the magnetic flux function, logio^ and the colour image shows = 2(1 + i<i)/c, where / is the electric current and 1^ 
is the displacement current. 



plasma is charge-separated and the force balance is achieved 
by means of a small parallel component of electric field (e.g. 
Polloway fc Price 198l"1 IMestel 1999^ . However the recent 
eclipse observations of the binary pulsar JO737-3039 allowed 
direct measurement of the particle density in the dead zone 
of one of the components | |Lyutikov fc Thompson 2005^ . It 
has turned out to be many orders of magnitude higher 
than the expected density of the charged-separated plasma. 
Whether such a high density is specific for binary systems 
only, as proposed in ( |Lyutikov fc Thom pson 2005), or typi- 
cal for single pulsars as well remains to be seen. In any case, 
the charge-separated dead-zones cannot be modelled within 
the MHD approximation where the required force balance 
can only be achieved by means of gas pressure. This would 
require the gas pressure to increase along the magnetic field 
line and reach a maximum in the equatorial plane. On the 
contrary, the magnetic pressure of dipolar magnetosphere 
decreases with distance as r -6 and thus the ratio of mag- 
netic to gas pressure has to decline even faster than r -6 . 
Given the requirement on the dead zone to remain magnet- 
ically dominated everywhere this implies very high pressure 
ratio near the star surface - much higher than our numer- 
ical scheme can accommodate. To overcome this problem 
one could consider the case where the star radius is only 2-3 
times smaller than the radius of its light cylinder but this 
would have a rather strong effect on the structure of the 
magnetosphere. This is why we have preferred a different 
strategy. 

Its idea is to reduce, somewhat arbitrarily, the dynam- 
ical effect of the centrifugal force on the motion of plasma 
near the star so that it does not destroy the nearly force-free 
equilibrium across the magnetic field lines within the dead 
zone. One way of achieving this is to reset the gas pressure 
and its rest mass density to some rather low 'target' values, 
p 3 and p s , within the dead zone every time step. At the same 



time one may reset the flow speed along the magnetic field 
lines to zero, just as it should be in equilibrium. However, 
the dead zone is not a well denned region in the case of 
time-dependent magnetospheres. Instead, one could apply 
the same procedure to a volume that is guaranteed to in- 
clude if not the whole dead zone then at least its inner part. 
For example, we used a sphere of radius r s < zoi c . This how- 
ever leads to another complication - in the open field region 
bounded by the sphere there emerges a strong rarefaction 
wave, so strong that the code crashes. Fortunately, this can 
be avoided if instead of resetting the target values one in- 
troduces a relaxation towards them with the relaxation time 
gradually increasing towards the boundary of the relaxation 
domain. We evolved p, p, and Vu according to the following 
equations 



dvn 



-6(2)(p-p«), 



= -&(2)(P-P«), 



(It 

dp 
It 

where 

= b (0).f(r), 
b(2){r,0) = 6(i) | cos f 
,., , , (r s -r)/(r B 

fin = ; 



if r < r 3 
if r > r 3 



(32) 
(33) 
(34) 

(35) 
(36) 

(37) 



where fo( j is constant, and r* is the star radius. As one 
can see, the relaxation time becomes infinite at r — r 3 . The 
results presented below correspond to r s = zui c . However, we 
have also tried smaller values of r a in order to verify that this 
does not lead to qualitatively different results (We discuss 
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Figure 3. Inner part of the solution at t = 55. Top left panel: The contours show the magnetic flux function, \P, and the colour image 
shows H^; Top right panel: The contours show the magnetic flux function, the arrows show the flow velocity, and the colour image 
shows the magnitude of the poloidal electric current density multiplied by r 2 ; Bottom left panel: The contours show the magnetic flux 
function, the arrows show the flow velocity, and the colour image shows the log 10 (wW 2 J B 2 ); Bottom right panel: The contours show 
the magnetic flux function, the arrows show the flow velocity, and the colour image shows B 2 . 



the effects of reducing r 3 in Section 5.) The dependence of 
6(i) on the azimuthal angle was introduce in order to reduce 
the possible adverse effect on the current sheet should it be 
formed inside the light cylinder. The actual value of b(o) is 
to be found by the method of trial and error. Finally, we use 
the following targets for the pressure and density 

p a = a (2 )p s c 2 , p s c = a (1) i3 2 , (38) 



axisymmetric domain (r, 9) = [0.1, 50] x [0, tt] and hence the 
star radius was set to r„ = 0.1. In order to speed up the 
calculation we started with a relatively low resolution grid, 
124 x 61, and then increased the resolution twice after the 
solution seemed to have reached a steady-state on the scale 
of several w% c . Hence the final grid had 496 x 244 cells. 



where am = 0.01 and a™ = 0.001. 

In these simulations, the computational grid covered the 
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4.3 Results 

The initial solution described a non-rotating magnetosphere 
with dipolar magnetic field, eq.JHJ, v* = 0, p — p s , and 
p — p s - At t = the star rotation is switched on and a tor- 
sional Alfven wave is emitted from its surface. As it propa- 
gates away larger and larger portion of the magnetosphere 
is set into rotation and develops an electric current system. 
This process is illustrated in figure [5] where the contours 
show the magnetic field lines and the colour image shows the 
distribution of = 2(1 + Id)/c, where I is the total elec- 
tric current and Id is the total displacement current through 
the circular contour of cylindrical radius zn = rsin#. Unity 
corresponds to = fiQ 2 c~ 2 . 

Behind the wave the solution gradually approaches a 
steady state. Figure [3] shows the inner region of this steady 
state solution at t — 55, the termination time of the simu- 
lations. The structure of magnetic field lines suggests that 
the dead zone extends all the way up to the light cylin- 
der. Beyond the light cylinder the poloidal magnetic field 
becomes radial, as expected |lngraham 1973| IMichel 197411 . 
Some magnetic field lines are closing up beyond the light 
cylinder but they do so within the equatorial current sheet 
due to finite artificial resistivity in the numerical scheme - 
it is easy to see the transition between the dead zone and 
the current sheet within which the magnetic field lines are 
highly stretched in the radial direction. The colour image in 
the top left panel of this figure shows the distribution of Hq. 
Since in a steady state the displacement current vanishes we 
have Htf, = 21 /c = A. This image confirms our conclusion 
on the extension of the dead zone which is seen in the im- 
age as a toroidal structure with H$ = (Indeed, within the 
dead zone the poloidal currents do not flow and H$ must 
be zero.) The jumps in occurring at the boundary of 
the dead zone and along the equator are indicators of thin 
sheets of return current. This image also shows, though not 
as clear, that H$ and hence I reach maximum amplitude 
at some finite distance from the equatorial plane. This in- 
dicates the presence of an additional layer of return current 
that surrounds the equatorial current sheet and the dead 
zone (Con topoulos et al. 1999| . 

The current sheets are most prominent in the right 
panel of figure [3] that shows the distribution the poloidal 
current density, J p , multiplied by r 2 . One can also see that 
the thickness of the current sheets near the Y-point is about 
0.5^0.6o7i c and that equatorial current sheet extends inside 
the light cylinder by approximately the same distance. 

From time to time it was suggested that particle iner- 
tia may actually become dynamically important near the 
light cylinder, thus rendering the force-free approximation 
as unsuitable. The argument develops like this. Provided 
the magnetospheric plasma simply co-rotated with the star 
its speed would exceed the speed of light beyond the light 
cylinder. This does not occur because of the rapidly increas- 
ing inertial mass of this plasma near the light cylinder. This 
leads to very strong centrifugal force that makes the plasma 
to flow across the light cylinder thus opening up the mag- 
netic field lines and sweeping them back. However, this not 
what occurs in our simulations. Indeed, the top left panel 
of figure |3] shows that, very much in agreement with the 
force-free models, is constant along the magnetic flux 
surfaces even when they cross the light cylinder. Moreover, 
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Figure 4. as a function of the magnetic flux function ^ on 
the sphere r = 1.1 at time t = 55. 

this conclusion is fully supported by the data presented in 
the bottom left panel of figure[3]that shows that the distribu- 
tion of wW 2 /B 2 , the quantity that can be used to describe 
the relative importance of the inertial effects as well as of 
the gas pressure. One can see that it remains very low ev- 
erywhere outside of the equatorial current sheet including 
the light cylinder. 

There is, however, one location in the force-free solution 
of Contopoulos et al.( where the inertial effects are ex- 

pected to become important. It is the so-called Y-point, that 
is the point where the dead zone approaches the light cylin- 
der. Indeed, since the dead zone co-rotates with the star 
then plasma particles attached to its field lines rotate with 
the speed reaching the speed of light at the y-point. More- 
over, one may argue that the divergence of magnetic field 
strength at the Y-point discovered by Gruzinov( 2005 1 also 
suggests a breakdown of the force-free approximation. This 
fully agrees with the data presented in the bottom right 
panel of figure [3] which shows the distribution of B 2 in our 
MHD solution. Although the plot shows a local minimum at 
(z, vj) — (0, 0.85) and then some growth of B 2 in the direc- 
tion of the Y-point this growth never develops. Moreover, in 
our solution the poloidal magnetic field lines approach the 
Y-point at an angle of ~ 50° to the equatorial plane (see 
figED> which is significantly lower than 77°. 3 predicted by 
Gruzinov( 2005 ) and quite close to the value of 56°. 5 found 
in Goodwin et al.(2004). In addition to the inertial effects, 
the relatively high thickness of the current sheets near the Y- 
point (the right top panel of figure^ and finite gas pressure 
also contribute to these discrepancies between our simula- 
tions and the force-free solution. Our results however do not 
prove that the asymptotic force- free solution for the Y-point 
is irrelevant. It seems possible that for sufficiently high mag- 
netization and small resistivity the exact MHD solution will 
first approach the force-free asymptote found by Gruzinov 
and then deviate from it on smaller scales. This, however, 
requires further investigation. 

Figure 0| shows H,/, as a function of Vl/ at r — 1.1 and 
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Figure 5. Left panel: The magnetic flux function, in the equatorial plane as a function of r at time t = 55. Right panel: The total 
energy flux through a sphere of radius r at time t = 55. 





Figure 6. The angular velocity distribution at t = 55 within the light cylinder at r = 0.2 (left panel) and at r ■ 
dot-dashed, dashed, and solid lines show the solutions with 61, 122, and 244 cells in the 6-direction. 



: 0.7 (right panel). The 



allows us to see these details of the current system more 
clearly. This distribution is very similar to the one given in 
figure 4 of Contopoulos et al. (1999). The minimum has 



1.02 



fiO, 



and A„ 



1.12 



The left panel of figure [K] shows the distribution of mag- 
netic flux function in the equatorial plane of the final solu- 
tion. As a consequence of the finite artificial resistivity in our 
numerical scheme the magnetic field lines continue to close 
up even beyond the line cylinder and this shows itself via a 
systematic decline of \t at r > 1. Thus, it is not so straight- 
forward to determine the fraction of opened field lines in this 
solution. One way to describe it quantitatively is by giving 
the value of the flux function exactly at (9 = n/2,r — 1). 
This gives us 



f yp ~ 1.37-^ 1.38^. 



(39) 



On the other hand, figure shows that the decline of $ 
significantly slows down at r > 2 where $ reaches the value 
of 



1.27 



(40) 



These numbers should be compared with the eSP ope „ / fj.il = 
1.23 in ContopouloslHIlEJ and Timokhin (|5nn5Jl . 1.27 in 
Gruzinov(2005), and 1.36 in Contopoulos et al.|l999). 

The right panel of figure shows the total flux of energy 
through the sphere of radius r that should be constant in 
a steady state solution. One can see that it is indeed more 
or less constant with the exception of 1 < r < 2 where the 
energy flux slightly increases. This increase is a permanent 
feature that arises because our scheme is not strictly conser- 
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Figure 7. The wind zone structure of a dipolar magnetosphcrc. 
The contours show the field lines of poloidal magnetic field. The 
colour image shows the distribution of log 10 W 



vative. In order to evaluate the spin-down power of the star 
we used the energy flux through the sphere of radius r = 1, 
thus ignoring the non-physical increase of total luminosity 
beyond r = 1. This gives us 



1.1 



(41) 



which is in a very good agreement with the result by Gruzi- 
nov(2005). 

Figure |S| shows the distribution of the angular velocity 
of magnetic field lines, Qf. In the exact steady-state force- 
free solution the magnetic field lines rotate with same angu- 
lar velocity as the star, 



Q f = Q. 



(42) 



However in our solution, Qf noticeably deviates from Q in- 
side the layer coincident with the current sheet between the 
open field lines and the dead zone (see fig|SJ ■ The most likely 
reason for this is the enhanced numerical resistivity in the 
current sheet. The comparison of solutions with different 
numerical resolution (fig|Sl shows that the thickness of the 
layer significantly decreases with resolution. However, the 
amplitude of the perturbations does not seem to depend on 
the resolution. 

Figure [7| shows that far away from the star the dis- 
tribution of poloidal field lines is very similar to the split- 
monopole one. This is exactly what was concluded in the pi- 
oneering papers by Ingrahamf 1973.) and Michel( 1974 J on the 
force-free magnetospheres. More recently, the centrifugally 
driven outflows in split-monopole magnetospheres have been 
studied within the cold-MHD approximation. According to 
these studies the Lorentz factor of centrifugally accelerated 
plasma at the fast critical point is 



W 



1/3 



(43) 



where a is the so-called magnetization parameter which is 



defined as the ratio of the Poynting flux density to the 
rest mass-energy flux density at the foot point of the mag- 
netic field line IBeskin 199711 . For outflows that are initially 
Poynting dominated a ~ a, where a is the ratio of the total 
energy flux density to the rest mass-energy flux density. In 
contrast to a, a can be measured at every point of the field 
line as it is constant along it. Figure|S|shows the distribution 
of 5", W, and the ingoing speed of the fast wave in the radial 
direction along the ray 9 — lrad. From these data we find 
that along this ray a ~ 77 and the fast point is located at 
r ~ 6.8 where W — 4.3. On the other hand, we can use the 
above value of a in order to calculate the Lorentz factor at 
the fast point directly from eg. 1431 - this give us W = 4.25. 
Thus, the split-monopole model provides quite a good model 
for the wind zone of aligned magnetic dipole. 

Another interesting feature of figure Q is the signifi- 
cantly higher value of the Lorentz factor of the outflow 
within the equatorial current sheet. This suggests that the 
mechanism of the flow acceleration is somewhat different 
there. The obvious suspect is heating due to resistive dissipa- 
tion of the electromagnetic energy. The high value of Lorentz 
factor in the current sheet suggests that its contribution into 
the global energy transfer can be quite significant. 

The left panel of figure [5] shows the evolution of the 
wind luminosities with the distance from the star by the 
end of simulations. The total luminosity is more or less con- 
stant apart from noticeable perturbations around r — 17, 
r — 30, and r — 45. The big bump around r — 45 is related 
to the leading front of the wind. The perturbations around 
r — 17 and r — 30 are related to the grid refinement events 
- each time the computational grid is refined the numeri- 
cal solution evolves to a slightly different state, most of all 
in the inner region of the computational domain. This trig- 
gers noticeable waves propagating away from the star. The 
electromagnetic luminosity, which is shown in figure |3 by 
the dashed line, gradually decreases with distance thus indi- 
cating the ongoing conversion of the electromagnetic energy 
into the hydrodynamic energy. This is supported by the evo- 
lution of the hydrodynamic luminosity of the wind which is 
shown by the dash-dotted line. One of the interesting fea- 
ture of this line is its rapid rise initial rise that suggests a 
particularly effective energy conversion. The nature of this 
conversion is clarified in the right panel of figure|3]that shows 
that the current sheet accounts for about 70% of the total 
hydrodynamic luminosity at r = 10 and this corresponds to 
about 15% of the total wind luminosity. This clearly points 
out at the Ohmic heating in the current sheet as the main 
source of the energy conversion. 



5 DISCUSSION 

There is not much to discuss in connection with our force- 
free simulations of pulsar magnetospheres. We simply have 
not been able to make the required progress using this frame- 
work because of the inability to handle the equatorial cur- 
rent sheet. For this reason we will focus in this section al- 
most entirely on the results of our MHD simulations and 
their possible implication for pulsar physics. 

One of the key goals of this study was to determine there 
the MHD equations allow stable, or quasi-stable, steady- 
state solutions for dipolar axisymmetric magnetospheres of 
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Figure 8. The variation of a number key parameters along the ray 8 = 1. Left panel: The local magnetization parameter ct; Middle 
panel: The wavespeed of the ingoing fast wave in the radial direction (the speed of light c=l); Right panel: The Lorentz factor of the 
wind. 
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Figure 9. Conversion of the electromagnetic luminosity into the hydrodynamic luminosity of the wind. Left panel: The total luminosity 
(solid line), the electromagnetic luminosity (dashed line), and the hydrodynamic luminosity (dot-dashed) of the wind at t = 55. Right 
panel: The angular distribution of the total luminosity at r = 10. The solid line and the dot-dashed line show the total luminosity and 
the hydrodynamic luminosity within the polar cone of angle 8 respectively; the dashed line shows the total flux density in the radial 
direction. 



neutron stars. This problem had become particularly inter- 
esting since the discovery a whole family of steady-state 
force-free solutions continuously parametrised by the loca- 
tion of their y-point IGoodwin et al. 20041 ITimokhin 2 005). 
Our results indicate the existence of a unique steady-state 
MHD solution to the problem and this solution is very close 
to the force-free stationary solution of the pulsar equation 
with the y-point located at the light cylinder - the orig- 
inal solution of pulsar found by Contopoulos et al.|l999). 
Why this solution is preferable to those whose y-point is 
located well inside the light cylinder has already been ex- 
plained in | |Contopoulos 2005^ . If resistivity is not vanish- 
ingly small then even if the initial solution had the dead 
zone buried well inside the light cylinder it would take only 
a finite time for the anti-parallel field lines of the current 
sheet between the y-point and the light cylinder to recon- 
nect and to form closed loops that become a part of the 
dead zone. The reconnection should also occur beyond the 



light cylinder but there the outflow is super- Alfvenic so the 
net outcome of the reconnection is likely to be the develop- 
ment of magnetic islands carried by the wind away from the 
star fUzdcnsky 2004). The rate of the reconnection depends 
on the actual resistivity in the current sheet, and in our 
simulations the resistivity is purely artificial. However, the 
ultimate outcome is unlikely to be different. Indeed, since 
the particle inertia on the closed field lines located well in- 
side the light cylinder of pulsar magnetospheres is extremely 
small there is no restoring force that would make this field 
lines to open up again. This is supported by the fact that 
the total electromagnetic energy of the force-free magneto- 
sphere is minimum when the Y-point is located on the light 
cylinder JTimokhin 20051 . 

At this point it makes sense to discuss whether the re- 
laxation procedure applied within r 3 = voi B (see Sec. 4. 2) 
could somehow promote the expansion of the dead zone to- 
wards the light cylinder. As we have already pointed out in 
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Figure 10. Comparison of solutions with relaxation domains of 
different sizes. The colour image shows H$ and the solid lines 
show the magnetic flux function for the model with r s = 0.7ro; c . 
The dashed lines show the magnetic flux function for the model 
with r s = zoi. 



Sec. 4.2 all relaxation times gradually increase towards in- 
finity as r — > r s . Thus, near the light cylinder the effect of 
the relaxation procedure on the flow is increasingly small. 
Moreover, the relaxation time for the gas pressure becomes 
infinite at 6 — tt/2. This would allow the build up of gas 
pressure required to support the equatorial current sheet 
should it exist within the light cylinder. However, in order 
to fully resolve this issue we carried out additional simula- 
tions with the radius of the relaxation sphere pushed down 
to r s = 0.7zui c . Figure [T01 allows to compare both solutions 
with regard to the location of the Y-point. The colour image 
in this figure shows the distribution of and the solid lines 
show the magnetic surfaces for the model with r 3 = 0.7ro; c 
whereas the dashed lines show the magnetic surfaces for the 
model with r s — zui c . As one can see more field lines open 
up for smaller r s and this is caused by the increased dy- 
namical role of the centrifugal force. However, the location 
of Y-point remains basically the same. Thus, our relaxation 
procedure cannot be considered as the reason for the dead 
zone extending all the way up to the light cylinder. 

Contopoulos( 2005 ) also pointed out that the open field 
lines of pulsar magnetospheres may rotate at a slower rate 
than the closed field lines of the dead zone due to the fi- 
nite potential gap of the polar cap. In particular, he spec- 
ulated about the possibility of a significant growth of the 
polar gap due to a sudden decrease of the particle injec- 
tion rate in the gap in order to explain the explosive phe- 
nomena like the December 27, 2004 burst in SGR 1806-20. 
Although, it is not clear what could cause such a sudden 
incident of "charge starvation" it was suggested long ago 
that a slow systematic evolution of the polar gap could re- 
sult from the gradual spin-down of the star ijSturrock 19711 
IRuderman fc S utherland 1975t . Thus, differentially rotating 



magnetospheres are not just interesting theoretical models 
but can be very much relevant for the real pulsars. 

Contopoulos( 2005 ) found force-free stationary numer- 
ical solutions for such magnetospheres in the simple case 
of a uniformly rotating polar cap and a uniformly rotating 
dead zone which was assumed to corotate with the star. (In 
fact, the dead zone may also have a potential gap that sep- 
arates electric charges of opposite sign but it is expected to 
be rather small ( |Holloway fc Price 1981^ .) These solutions 
have two light cylinders - a smaller one for the dead zone 
and and a larger one for the open field lines and the dead 
zone extends all the way up to its light cylinder. Contopou- 
los(2005) argued that although there existed solutions with 
a smaller dead zone they were not sustainable due to recon- 
nection in the part of the equatorial current sheet that runs 
between the y-point and the light cylinder of the dead zone. 

In fact, the numerical models constructed by Contopou- 
los(2005) are likely to be globally unstable to reconnection 
too. Indeed, in these models the equatorial current sheet 
continues inside the light cylinder of open magnetosphere 
and magnetic reconnection occurring in this region should 
lead to creation of new closed field lines. However, this case 
is somewhat more involved. Let us imagine that such re- 
connection has indeed occurred. The field lines that have 
just closed down are now beginning to spin-up and coro- 
tate with the dead zone. However, they extend beyond the 
light cylinder of the dead zone and for this reason they can 
not corotate with it - as they spin-up they begin opening 
up again. Once have been opened up they begin to slow 
down, thus creating conditions for the next closing down 
event, and so on. This simple analysis suggests that such 
differentially rotating magnetospheres cannot be stationary 
and have to develop oscillations. The typical time scale of 
such magnetospheric oscillations seems to be determined by 
the spin-up time that must be comparable with the time 
required for an Alfven wave to cross the distance between 
the light cylinder of the dead zone and the star back and 
forward. Since the Alfven speed is relativistic and the mag- 
netic field has a significant radial component the crossing 
time has to be comparable with with rotational period of 
the star. The reconnection rate is more likely to effect the 
amplitude of these oscillations rather than their time-scale, 
a quicker reconnection leading to a larger fraction of mag- 
netic field lines involved in this process of closing down and 
opening up. We wont to speculate that these oscillations 
may be relevant to the origin of the sub-pulses of radio pul- 
sars, e.g. ((Manchester & Tay lor 19770 . For periodic magne- 
tospheric oscillations we would have a phenomena reminis- 
cent of beating waves - this may explain the so-called drift- 
ing sub-pulses. However, one may also expect quasi-periodic 
and even chaotic oscillations and they would results in a 
much more complicated behaviour of sub-pulses. 

Our MHD solution has a number of interesting features 
that could not possibly be found in the ideal force-free so- 
lution of Contopoulos et al.(1999). Some of them do not 
depend much on the details of resistivity like, for example, 
the centrifugal acceleration of the wind outside of the cur- 
rent sheet. Others, like the Ohmic dissipation and the wind 
acceleration in the equatorial current sheet, do and we have 
to exercise a reasonable degree of caution when interpret- 
ing them - MHD simulations with only artificial resistiv- 
ity can provide at most qualitatively correct description of 
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such features. However, it is interesting that the dissipation 
of electromagnetic energy in the equatorial current sheet 
has already been considered as a promising explanation of 
the high luminosity gamma-ray emission from young pul- 
sars l |Lyubarskii 1996| IKirk et al. 20 02 ) . Given the fact that 
these gamma-rays carry away a significant fraction of the 
spin-down power, up to 10% in the most extreme exam- 
ples ({Thompso n 200 1) , and that pulsar magnetospheres are 
highly magnetically dominated, an efficient dissipation of 
Poynting flux somewhere in the magnetosphere is needed to 
explain these observations, and the equatorial current sheet 
is one of the most natural locations for such dissipation. 
In particular, Lyubarskii (1996) argued that this high en- 
ergy emission originates from the current sheet just beyond 
the dead zone (This is exactly where our simulations show 
most effective Ohmic dissipation.) The polarizational obser- 
vations of the optical emission from the Crab pulsar support 
this idea. This emission is found to be polarised parallel to 
the rotation axis at the peaks of the pulses and between 
the pulses (within the pulses the polarisation vector rotates, 
which may be attributed to the rotation of the magneto- 
sphere). This shows that the magnetic field of the emitting 
region is predominantly perpendicular to the rotation axis 
as it should be if this radiation is generated in the equatorial 
current sheet. 
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